%------------------------------------------------------------------------%
%              Comunicaciones digitales y analógicas                     %
%                      Trabajo Práctico N 1                              %
%                                                                        %
%                       Eric Lionel Koplin                               %
%------------------------------------------------------------------------%
close all;clear all;clc;
%% Probabilidad de error de una M-PAM
%---------------------------------------------------------------Parámetros
% Parametros de ajuste de la simulacion
%

Pes=zeros(4,13);
Peb=zeros(4,13);

for modindex=1:4    % distintos moduladores
    for i=1:4       % distintas cantidades de símbolos en el modulador

        M = 2^i;        % cantidad de simbolos
        k = log2(M);    % bits por simbolo
        SNRbdB = 2:14;  % SNR por bit en dB
        nsimbs = 2e4 ;  % cantidad de simbolos que para cada valor de SNRbdB en el Montercarlo

        %----------------------------------- Conversiones e inicializaciones varias
        %
        SNRbveces = 10.^(SNRbdB/10); % SNR por bit en veces
        SNRdB = SNRbdB + 10*log10(k);% SNR en dB
        SNR = 10.^(SNRdB/10);% SNR en veces

        Eavg = 1 ;
        Pes = zeros(size(SNRbdB));% probabilidad de error por simbolo simulada
        Peb = zeros(size(SNRbdB));% probabilidad de error de bit simulada

        %--------------------------------------------------------P(error) teorica :
        %
        % Probabilidad de error de bit teorica de una 4−PAM
        % Uso un puntero a una funcion cuyo parametro es SNR (en veces)
        Pet = @(SNR) 2*(1-1/M) * qfunc(sqrt(3*SNR/(M^2-1))) ;

        %-----------------------------------------------------Simulacion Montecarlo
        
        switch modindex
            case 1
            %--------------------PAM
            % modulador
            modulator = comm.PAMModulator(M,'SymbolMapping','Gray',...
                                      'NormalizationMethod', 'Average Power',...
                                      'AveragePower',Eavg);
            %demodulador
            demodulator = comm.PAMDemodulator(M,'SymbolMapping','Gray',...
                                          'NormalizationMethod' , 'Average Power' , ...
                                          'AveragePower',Eavg);
            case 2
            %--------------------PSK
            %Recordar la energia está normalizada, Es=Eavg=1
            % modulador
            modulator = comm.PSKModulator(M,'SymbolMapping','Gray');% Eavg=1
            %demodulador
            demodulator = comm.PSKDemodulator(M,'SymbolMapping','Gray');
            
            case 3
            %--------------------QAM
            % modulador
            modulator = comm.RectangularQAMModulator(M,'SymbolMapping','Gray',...
                                      'NormalizationMethod', 'Average Power',...
                                      'AveragePower',Eavg);
            %demodulador
            demodulator = comm.RectangularQAMDemodulator(M,'SymbolMapping','Gray',...
                                          'NormalizationMethod' , 'Average Power' , ...
                                          'AveragePower',Eavg); 
            
            %case 4
            %--------------------FSK
            %Recordar la energia está normalizada, Es=Eavg=1
            % modulador
            %modulator = comm.FSKModulator('ModulationOrder',M,'SymbolMapping','Gray',...
            %                            'SamplesPerSymbol',1);% Eavg=1
            %demodulador
            %demodulator = comm.FSKDemodulator('ModulationOrder',M,'SymbolMapping','Gray',...
            %                            'SamplesPerSymbol',1);                          
        end
        %--------------------AWGN Channel
        %hAWGN = comm.AWGNChannel('NoiseMethod','Signal to noise ratio (Eb/No)',...
        %    'EbNo',SNRbdB(1),'SignalPower',Eavg,'BitsPerSymbol',k);
        %hAWGN = comm.AWGNChannel('NoiseMethod', ...
        %                         'Signal to noise ratio (SNR)','SNR',0);
        %--------------------Create an error rate calculator
        %hError = comm.ErrorRate;

        %recordar hacer sigma/2 en el teorico tambien
        simbtx = randi(M,nsimbs*length(SNRbdB),1)-1;% simbolos a transmitir (en 1:M);
        
        if modindex~=4
            simbmod = step(modulator,simbtx); % simbolos modulados 
            %simbmod_ruido = step(hAWGN,reshape(simbmod,nsimbs,length(SNRbdB)));%en cada columba una canal distinto
            sigma = sqrt(Eavg./SNR); % desvio de ruido
            ruido = complex(randn(nsimbs,1),randn(nsimbs,1))*sigma;
            simbmod_ruido = simbmod + reshape(ruido,nsimbs*length(SNRbdB),1);

            simbrx= step(demodulator,reshape(simbmod_ruido,nsimbs*length(SNRbdB),1)) ; % simbolos recibidos
            %errores = step(hError, simbtx, simbrx);
        else %FSK -- modulacion ortogonal
            simbmod = zeros(length(simbtx)*M,1); % en grupos de M las coordenadas m-dimensionales
            simbmod([0:nsimbs*length(SNRbdB)-1]'*M+simbtx+1)=1;
            %
            sigma = sqrt(Eavg./SNR); % desvio de ruido
            ruido = complex(randn(nsimbs*M,1),randn(nsimbs*M,1))*sigma;% nM x length(sigma)
            
            simbmod_ruido = reshape(simbmod,nsimbs*M,length(sigma)) + ruido;

            [modulos simbrx]=max(abs(eye(M)*reshape(simbmod_ruido,M,[]))); % simbolos recibidos
            simbrx = simbrx'-1;
        end
        %cada columna es un canal con distinta snr y cada elemento es una señal

        Pes(i,:) = ones(1,nsimbs)*reshape(simbtx ~= simbrx,nsimbs,length(SNRbdB))/nsimbs ;

        % Ahora analizo los errores de bit (el Gray lo hizo el modulador )
        %
        bitstx = de2bi(simbtx,k) ; % convierto los simbolos transmitidos a bits
        bitsrx = de2bi(simbrx,k) ; % convierto los simbolos recibidos a bits
        % comparo los bits y cuento cuantos difieren
        Peb(i,:) = ones(1,nsimbs)*reshape(xor(bitstx,bitsrx)*ones(k,1),nsimbs,length(SNRbdB)) /k/nsimbs;

        % Graficos
        %
        % Probabidad de error de simbolo
        figure(2*modindex-1)
        semilogy(SNRbdB,Pet(SNR),'-dr','linewidth',1.5);hold on
        semilogy(SNRbdB,Pes(i,:),'-o','linewidth',1.5);
        grid on
        box on
        xlabel('SNR por bit (dB)') ;
        ylabel('Pe')
        legend(' Curva teorica','Montecarlo')

        % Probabilidad de error por bit
        figure(2*modindex)
        semilogy(SNRbdB,Pet(SNR)/k,'-dr');hold on
        semilogy(SNRbdB,Peb(i,:),'-o','linewidth',1.5);
        grid on
        box on
        xlabel('SNR por bit (dB)') ;
        ylabel('Pb')
        legend('Curva teorica','Montecarlo')
    end
end